SUBROUTINE ARRGATHERWRI_MPI(array,i0,i1,j0,j1,k0,k1, &
& i0_out,i1_out,j0_out,j1_out,k0_out,k1_out,nx1,ny1,ns1, &
& varname)

! Subroutine ARRGATHERWRI_MPI gathers the portions of the real type array from
! all MPI processes at one process and outputs the whole array in ASCII file.
! Tha array is assumed to be decomposed in two horisontal dimensions.

!use MPI
use MPI_VARIABLES, only : &
& arrgatherwri_indic, &
& nxpi, nxpe, nypi, nype, &
& comm3d, tag, size_MPI, & 
& rank_comm3d, ierr_MPI, &
& parallel
use TIMING_MOD, only : &
& proc_time

implicit none

include 'mpif.h'

! Input variables
real(8), intent(in) :: array(i0:i1,j0:j1,k0:k1)

integer(4), intent(in) :: i0, i1
integer(4), intent(in) :: j0, j1
integer(4), intent(in) :: k0, k1 ! Must be the same for all processes
integer(4), intent(in) :: i0_out, i1_out
integer(4), intent(in) :: j0_out, j1_out
integer(4), intent(in) :: k0_out, k1_out
integer(4), intent(in) :: nx1, ny1, ns1

character(len=*), intent(in) :: varname

! Local variables

real(8), allocatable :: global_array(:,:,:)
real(8), allocatable :: buf(:), bufrecv(:,:)
real(8) :: pr_time, pro_time

integer(4) :: i, j, k, ijk, ijk2, ijk_max, irank
integer(4) :: tag1 = 1
integer(4) :: nunit = 345
integer(4) :: status_recv(1:MPI_STATUS_SIZE,1:size_MPI-1)
integer(4) :: status_send(1:MPI_STATUS_SIZE)
integer(4) :: request_recv(1:size_MPI-1) 
integer(4) :: request_send

character(len=5) :: nproc_char
character(len=100) :: format_char
character(len=100) :: filename

logical, allocatable :: ok_global_array(:,:,:)

!pr_time = MPI_WTIME()
call cpu_time(pr_time)

if (parallel) then
! Sending the portion of the process to the process with rank 0        
  if (rank_comm3d /= 0) then
    ijk2 = (nxpe(rank_comm3d)-nxpi(rank_comm3d)+1) * &
    &      (nype(rank_comm3d)-nypi(rank_comm3d)+1) * (k1 - k0 + 1)
    allocate (buf(1:ijk2))
!    print*, 'Sended-1'

    ijk = 0
    do k = k0, k1
      do j = nypi(rank_comm3d), nype(rank_comm3d)
        do i = nxpi(rank_comm3d), nxpe(rank_comm3d)
          ijk = ijk + 1
          buf(ijk) = array(i,j,k)
        enddo
      enddo
    enddo
!    print*, 'Sended-2'

    tag = 0
    call MPI_ISEND(buf,ijk2,MPI_DOUBLE_PRECISION, &
    & 0,tag,comm3d,request_send,ierr_MPI)

    call MPI_WAIT(request_send, status_send, ierr_MPI)
    deallocate (buf)
!    print*, 'Sended-3'
  else
! Receiving the portions of the array from all the processes 
! at the process with rank 0
    allocate (global_array(0:nx1,0:ny1,k0:k1))
    allocate (ok_global_array(0:nx1,0:ny1,k0:k1))
    ok_global_array(:,:,:) = .false.
    ijk_max = 1
    do irank = 1, size_MPI-1
      ijk_max = max(ijk_max, &
      & (nxpe(irank)-nxpi(irank)+1)*(nype(irank)-nypi(irank)+1)*(k1-k0+1))
    enddo
    allocate (bufrecv(1:ijk_max,1:size_MPI-1))

    do irank = 1, size_MPI-1
      ijk = (nxpe(irank)-nxpi(irank)+1)*(nype(irank)-nypi(irank)+1)*(k1-k0+1)
      call MPI_IRECV(bufrecv(1,irank),ijk,MPI_DOUBLE_PRECISION, &
      & irank,tag,comm3d,request_recv(irank),ierr_MPI)
    enddo

    do irank = 1, size_MPI-1
      call MPI_WAIT(request_recv(irank), status_recv(1,irank), ierr_MPI)
      ijk = 0
      do k = k0, k1
        do j = nypi(irank), nype(irank)
          do i = nxpi(irank), nxpe(irank)
            ijk = ijk + 1
            global_array(i,j,k) = bufrecv(ijk,irank)
            ok_global_array(i,j,k) = .true.
          enddo
        enddo
      enddo
    enddo

    deallocate (bufrecv)

    do k = k0, k1
      do j = nypi(0), nype(0)
        do i = nxpi(0), nxpe(0)
          global_array(i,j,k) = array(i,j,k)
          ok_global_array(i,j,k) = .true.
        enddo
      enddo
    enddo

  endif

else

  allocate (global_array(0:nx1,0:ny1,k0:k1))
  allocate (ok_global_array(0:nx1,0:ny1,k0:k1))
  global_array = array(0:nx1,0:ny1,k0:k1)
  ok_global_array(:,:,:) = .true.

endif

if (rank_comm3d == 0) then

  do k = k0, k1
    do j = 0, ny1
      do i = 0, nx1
        if (.not.ok_global_array(i,j,k)) then
          write(*,*) 'Error in ARRGATHERWRI: STOP'
          STOP
        endif
      enddo
    enddo
  enddo

  call CHECK_UNIT0(nunit)
  write(nproc_char,'(i5.5)') size_MPI
  filename = 'array_3d_'//nproc_char//varname
  open (unit=nunit,file=filename(1:len_trim(filename)), &
  & status='unknown')

!  format_char = '(3i6,f30.20)'
  format_char = '(3i6,E40.30)'
  do k = k0_out, k1_out
    do j = j0_out, j1_out
      do i = i0_out, i1_out
        write (nunit,format_char) i, j, k, global_array(i,j,k)
      enddo
    enddo
  enddo
  close (nunit)
  deallocate (global_array)

endif

!call MPI_BARRIER(comm3d,ierr_MPI)

!proc_time(arrgatherwri_indic,2) = proc_time(arrgatherwri_indic,2) + &
!& MPI_WTIME() - pr_time
call cpu_time(pro_time)
proc_time(arrgatherwri_indic,2) = proc_time(arrgatherwri_indic,2) + &
& pro_time - pr_time


END SUBROUTINE ARRGATHERWRI_MPI

